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CN ■ A computer simulation model is used to study the density profile 

and flow of a miscible gaseous fluid mixture consisting of differing con- 
stituent masses (tiia = wib/3) through an open matrix. The density 
profile is found to decay with the height oc exp(— m^mh), consistent 
with the barometric height law. The flux density shows a power-law 
increase oc (p c — p)^ with [i ~ 2.3 at the porosity 1 — p above the pore 
'"^3 ' percolation threshold 1 — p c . 
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1 Introduction 



Understanding the flow of a complex gaseous mixtures, sedimentation, and 
evolution of density profiles of its constituents in geo-marine systems and 
near-surface ecological environments is becoming increasingly important Jl|, 
|J. There are a number of examples: (i) High density brines associated with 
salt tectonics in large salt provinces (e.g. the Gulf of Mexico) have been 
breaching the seafloor and forming pools toxic to native flora and fauna Q. 
This involves flow of a mixture of fluids with different densities due to dif- 
ferent salt content and temperature, (ii) A mixture of air and radon flux 
through unsaturated soils within the upper few meters of the land surface. 
Radon 222, one of the intermediate products of the decay of uranium 238 to 
lead 206, is an odorless, radioactive gas (with half life 3.8 days), and is com- 
mon in many soils and rocks. Because radon is about 8 times more dense than 
air, and is relatively inert Q , it easily penetrates porous building materials in 
ground floors and basements, especially when pressure gradients are created 
by central heating systems 0. The U.S. Environmental Protection Agency 
(EPA) estimated in 1986 that 5,000 to 20,000 persons in the United States die 
of lung cancer each year from inhaling radioactive radon decay products in 
homes and buildings [H. (in) Evidence of methane hydrate formation below 
the ocean floor and in sub-ocean bottom in mud-volcano involves the flow 
and sedimentation of complex gas and fluid mixtures with dissimilar masses 
H [7], || . Studies of flow and density distribution of a mixture of miscible gas 
and fluids through porous media are therefore highly desirable. 

Systematic studies based on the field measurements of flow and density 
profile of gas and fluid constituents in geomarine environment are severely 
limited U due to uncontrollable changes in geothermal parameters and mor- 
phological variations. Thus, computer simulation studies remain one of the 
most viable tools to probe such difficult issues as flow H, [TO, ITTJ and density 



profile |TJ|. Incorporating all details (multiple constituents and their charac- 
teristics) even with a coarse grained host matrix i.e. open porous media with 
appropriate concentration gradient and pore distributions is a challenging 
issue. Lattice gas and particulate methods in general (Boltzmann, Cellular 



Automata, Ising (interacting), etc.) [|T^, [T^j have been used in diverse ap- 
plications of fluid flow. In study of the density profile and fluid flow of an 
interacting fluid mixture through porous media, a direct application of tra- 
ditional hydrodynamics approaches fl5| , |16|1 becomes intractable; apart from 
a major problem of enormous boundary conditions in such a porous medium 
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(percolating system), it is not clear how to include interaction or reaction be- 
tween the fluid components in hydrodynamic equations. Interacting lattice 
gas |12| may be a simple approach to initiate probing such difficult issues. 

Very recently we studied the flow of a fluid described by particles, say 
of type "A" through an open porous medium using a computer simulation 
model |T7[ . The porous medium is generated by a random distribution of 
sediment particles of concentration p on a simple cubic lattice. The bottom 
layer of the matrix is connected to a source of mobile fluid particles ("A"). 
As soon as a bottom site becomes empty, it is immediately filled by a particle 
"A" from the source. Particles in the bottom layer are not allowed to move 
below this plane due to presence of abundant source particles. On the other 
hand, the particles can escape the system from the top if they attempt to 
move to the higher layer. In this concentration driven system, the flux density 
shows a power-law decay with the porosity near percolation threshold. The 
steady-state density profile of fluid particles depends systematically on the 
barrier concentration p. 



In this article, we extend our previous studies [ff7| to a miscible two- 
component system consisting of constituents, say "A" and "B" with dissimi- 
lar masses (m a and nib, rrib = 3m a ). The model is described in the following 
section 2. We incorporate the effect of gravitational potential energy in mov- 
ing the particles and allowing them to escape the system from the bottom 
layer as well. The injection probability of particles A and B at the bottom 
remains equal. The results are presented in section 3 with a conclusion in 
section 4. 



2 Model 

Each site on a simple-cubic Lx L x L lattice, with L up to 250, can be in one 
of four states: occupied by an A particle, occupied by a B particle, empty 
(0), or a barrier site. Nearest-neighbor particles interact with energy J such 
that A and B mix well: J(A, A) = J(B, B) = -J(A, B) = -J(B, A) = 
— J(A,0) = —J(B,0) = 1 where negative J means attraction and positive 
J means repulsion. The immobile barriers exert no force on the particles 
except to prevent them from occupying the barrier site. The energy thus is 

i k 
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where % runs over all particles, k over all neighbor sites of i, and I and K are 
the corresponding site variables (A, B, barrier, or 0). 



Time development for one 250 fl 3 lattice; t-10 (+), 100 (x), 1000 (*), 5000 (full sq.), 10000 (open sq.) 
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Time development for one 250 fl 3 lattice; t=10 (+), 100 (x), 1000 (*), 5000 (full sq.), 10000 (open sq.) 
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Figure 1: Crossover from constant to exponential density profiles, with 
increasing times t as given in headline, for A-particles (part a) and B-particles 
(part b). No sites are barriers in this simulation. 

A and B particles can move with a Metropolis probability exp(—AE/ksT) 
to neighboring sites; here AE is the energy change associated with this move 
and ksT = 5 is the thermal energy. If a particle on the lowest plane moves 
upwards or horizontally, we inject a new particle A or B (with equal proba- 
bility) onto the vacated site. If a particle on the lowest plane wants to move 
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downward, it drops out. A particle moving upward from the highest plane 
(relatively a rare event) is lost, without any injection from the top. Periodic 
boundary conditions are applied in the horizontal directions. In this way, 
an oo x oo x L plate is approximated, with new material injected from the 
bottom at a nearly constant rate. 

Gravity pulls down the particles through an energy m/c^T where the lat- 
tice constant and the gravitational constant are incorporated into the dimen- 
sionless mass m = 0.1 for A and m = 0.3 for B particles. Thus the barometric 
height law gives an equilibrium density oc exp(— mh) as a function of height 
h, 1 < h < L. 

One time step is an attempt to move each particle once (on average) 
through random sequential updating; it does not matter much if instead we 
enforce exactly one attempt per particle for each time step. For L = 30 
we used up to t — 10 6 time steps, without seeing any long-time effects; for 
larger L (50 to 250) typically t = 10 4 to 10 5 gave equilibrium. Initially, the 
lattice is occupied homogeneously with a low concentration of particles, half 
A and half B. Our computer program allows many more choices for inter- 
actions and boundary effects and is developed to investigate many different 
systems (details are available from ras.pandey@usm.edu). One diffusion at- 
tempt, without barrier sites, took about 0.5 microseconds on one Cray-T3E 
processor. 

3 Results 

Fig.l shows for our largest lattices how the initial constant density profiles 
change with time, starting from the two boundaries, into a roughly exponen- 
tial decay for intermediate heights. The equilibrium density profiles, Fig.2, at 
intermediate densities are consistent with the barometric height law shown as 
straight lines in these semilogarithmic plots, both with and without barrier 
sites. The nearest-neighbor correlation functions (not shown), i.e. the num- 
ber of A or B particles surrounding a particle at height h, decay qualitatively 
similar to the density profiles. 

Fig. 3 shows as a function of the barrier concentration p the system's 
permeability, defined as the net flow (per unit time and unit cross-sectional 
area) at the top or bottom surfaces divided by the injection rate at the lowest 
plane. The straight line in this log-log plot suggests a power law similar to 
that of the electrical conductivity in random resistor networks. Fig.4 shows 
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Density of A (+) and B (x) particles, 64 lattices 100 A 3, t = 10 A 5, barrier concentration 1/2 
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Density of A(+) and B(x) particles, 64 lattices 100 A 3, t = 20000, no barriers 
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Figure 2: Equilibrium density profiles with (part a) and without (part b) 
half the sites occupied by barriers. The straight lines are oc exp(— mh). 



the density profiles of the B particles for the same simulations. We see there 
an exponential decay, followed by a plateau whose value strongly increases 
with increasing p, suggesting a trapping of particles by the barrier sites. We 
have also looked at the velocity and instantaneous velocity profiles of the two 
constituents which are consistent with our expectations, i.e., more mobility 
toward the top. 
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Figure 3: Double-logarithmic plot of permeability versus p c — p where p c = 
0.6884 is the percolation threshold and p the barrier concentration. For 
p > p c there is no infinite connected set of fluid sites between the barriers 
and thus the permeability vanishes. The straight line has a slope 2.3. 

4 Conclusion 

A computer simulation model is used to study the flow rate, sedimentation, 
and density profile of a mixture of miscible particle systems for a range of 
porosities above the pore percolation threshold in an open porous medium. 
In our concentration driven system, the steady-state density profiles for both 
A and B fluid are reached. Both density profiles show the well-known expo- 
nential decay cx exp(— mA(B)h) with height h. The effect of mass difference 
is vividly illustrated in the density profiles: while for our lattice sizes the 
density of A particles continues to decay up to the top plane the density 
of B particles already saturates at some low value. The saturated density 
of B particles increases on decreasing the porosity 1 — p - we speculate this 
saturation is due to trapping of particles in the pores. 
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Barrier concentrations 0.0 (+), 0.1 (x), 0.2 (*), 0.3 (open sq.), 0.4 (full sq.), 0.5 (open o), 0.6 (full o) 
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Figure 4: Equilibrium density profiles of B-particles for various barrier 
concentrations (same simulations as for Fig. 3). 

The flux of particles A at the bottom becomes equal to the outward flux 
from the top in steady state. The flux density of particles A decays with 
porosity with a power-law oc (p c — p)^ with /i ~ 2. 
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